Wave optics of imaging with contact ball lenses

Recent progress in microspherical superlens nanoscopy raises a fundamental question about the transition from super-resolution properties of mesoscale microspheres, which can provide a subwavelength resolution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim \lambda /7$$\end{document}∼λ/7, to macroscale ball lenses, for which the imaging quality degrades because of aberrations. To address this question, this work develops a theory describing the imaging by contact ball lenses with diameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$30<D/\lambda <4000$$\end{document}30<D/λ<4000 covering this transition range and for a broad range of refractive indices \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.3<n<2.1$$\end{document}1.3<n<2.1. Starting from geometrical optics we subsequently proceed to an exact numerical solution of the Maxwell equations explaining virtual and real image formation as well as magnification M and resolution near the critical index \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\approx 2$$\end{document}n≈2 which is of interest for applications demanding the highest M such as cellphone microscopy. The wave effects manifest themselves in a strong dependence of the image plane position and magnification on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D/\lambda $$\end{document}D/λ, for which a simple analytical formula is derived. It is demonstrated that a subwavelength resolution is achievable at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D/\lambda \lesssim 1400$$\end{document}D/λ≲1400. The theory explains the results of experimental contact-ball imaging. The understanding of the physical mechanisms of image formation revealed in this study creates a basis for developing applications of contact ball lenses in cellphone-based microscopy.

near-fields under coupling with whispering gallery modes (WGMs) in microspheres 27 , and extreme curvature of the wavelength-scale microspheres 23 are being considered to explain experimentally observed super-resolution values.
The purpose of this work can be seen as filling the gap between the classical imaging approach based on geometrical optics with aberrations included in analysis for sufficiently large ball lenses on the one hand and a limiting case of contact, wavelength-scale microspheres for which introducing aberrations is conceptually difficult on the other hand. We focus specifically on objects near the optical axis where the spherical aberration and defocus aberration dominate. We consider the imaging by ball lenses with the index of refraction 1.3 < n < 2.1 and radius 15 < R/ < 2000 which are placed in contact with the objects. Such imaging can be viewed as an extension of MSI towards using submillimeter-and millimeter-scale ball lenses which are much more convenient for developing portable and lightweight imaging systems not requiring bulky and expensive microscope objectives with heavy microscope stands. In addition, millimeter-scale ball lenses provide larger field-of-view compared to microspheres. The perfect example of such a system is a cellphone operating in combination with a contact ball lens 31 illustrated in Fig. 1a. The interest to similar systems operating with noncontact ball lenses emerged in the last decade [32][33][34][35][36][37] . However, some problems such as insufficient magnification, spherical aberration, and pincushion distortions were found to be limiting factors for developing this technology with the best resolution about 1.5 µ m reported up to date.
In our theoretical analysis we make a gradual transition from the paraxial to nonparaxial ray tracing and, finally, to wave optics (WO) based on a numerical solution of the Maxwell equations which provide an accurate description of imaging by the contact ball lenses with various refractive indices and sizes. As expected, our WO results converge to that of GO for larger ball lenses, while they also show a large number of novel properties taking place in the intermediate range of sizes. Using a two-dimensional (2-D) model, we study the effect of focal shift for contact ball lenses for the first time. It is shown that it plays a significantly larger role for real imaging taking place for 1.9 < n < 2.1 in the size range 50 < R/ < 2000 , compared to the virtual imaging observed for 1.3 < n < 1.9 . We also show that the effect of focal shift becomes more significant for smaller ball lenses leading to a dramatic reduction of the image magnification compared to the GO predictions. It is shown that the resolution improves for smaller ball lenses reaching deeply subwavelength value ∼ /2 at R/ 100 in the vicinity of n ≈ 2 . The index variation about n ≈ 2 affects greatly the magnification but not the resolution. Using diffraction integrals we clarify the role of defocus and spherical aberration on the shift of the image plane and also derive a simple formula describing the focal shift in a wide range of parameters , n, and R in a good agreement with the exact numerical solution of the Maxwell equations. We show that the WO results for the image plane shift agree well with the experimental cellphone microscopy of a Siemens star through a contact R = 1 mm ball lens at three different values of its refractive index. Finally, we experimentally demonstrate that cellphone microscopy based on using ball lenses with n sufficiently close to the critical index n ≈ 2 allows reaching ∼ 0.9 µ m resolution values in a resonable agreement with theoretical predictions.

Geometrical optics of ball lenses
Geometry and typical experimental situation. The use of ball lenses with the index of refraction n ≈ 2 can provide very large magnification 31 . The experimental setup for imaging through a ball lens is shown in Fig. 1a. The measurement details are described in Methods. The imaged object was a Siemens star and the real images were recorded by a cellphone camera, see Fig. 1b. The variation of the refractive index near n ≈ 2 was obtained due to the material dispersion of the glass and the selection of the operating wavelength using optical filters. The measured magnification can be very high, in the range 20-50, but depends strongly on the refractive Scanning electron microscopy (SEM) image of a Siemens star and its cellphone camera images through the ball (LASFN35 glass, R = 1 mm) using optical filters at the specified wavelengths. (c) Geometrical optics of imaging with a ball lens for an object at the surface. A virtual image is created for n r < 2 and a real image for n r > 2 . The black arrows indicate the displacement of the image as n r → 2.
Paraxial geometrical optics. Let us now analyze the experimental results on the ball-lens imaging within the framework of paraxial GO. The application of GO is rather natural here since the size of the ball is much larger than the wavelength. We consider a ball with radius R and refractive index n 2 surrounded by medium with n 1 , see Fig. 1c. An object is located next to the surface of the ball so that the gap is significantly smaller than the wavelength. The variation of the gap g(y) due to curvature can be estimated as g(y) = y 2 /(2R) , where y is the displacement from the optical axis. Assuming a subwavelength condition g < /2 gives the maximum object size L = 2 √ R for which the subwavelength gap approximation holds. Taking, for example, R = 1 mm and = 0.5 µ m yields L ≈ 40 µ m. Thus, the subwavelength gap approximation can easily hold for various samples of interest, such as used in the experiments, see Fig. 1b. Thus, the first interface does not participate in the image formation. The emission diagram formed inside the lens is transformed by the refraction at the second interface. The curvature of the first interface, however, can play a role in increasing the resolution due to the outcoupling of evanescent fields but this is likely to take place for lenses with sizes in the micron range. Ray 1, which originates from a point at distance h ≪ R from the optical x axis and propagates parallel to it, upon refraction bends at an angle γ 1 = (n r − 1)h/R , where n r = n 2 /n 1 , while ray 2, which passes through the center of the ball, has γ 2 = h/R . The equations for the two refracted rays allow finding their intersection x i and resulting object magnification M: This equation was also used to evaluate magnification in virtual imaging 10,24 and later generalized to the case of a finite gap between the object and the ball lens 12 . The intersection is virtual ( x i < 0 ) if γ 1 < γ 2 or n r < 2 and real ( x i > 0 ) in the opposite case. According to Eq. (2), the distance |x i | and magnification M increase rapidly as n r → 2 . Since M ∝ x i , we can limit ourselves to finding the image location x i for a point source on the optical axis.
The image plane positions x i /R predicted by the paraxial GO for the experimental parameters using Eq. (2) are significantly larger than the measured values, see Table 1. Moreover, the difference increases drastically as n r approaches 2. In general, the paraxial approximation is applicable only if the image is formed by rays in a small angular range, often limited by apertures in practice. The failure of the paraxial GO can be explained by the fact that the imaging here was performed using the full angular range offered by the ball lens.
Ray tracing beyond the paraxial approximation. Let us now go beyond the paraxial approximation and consider various rays originating from a point source located at x = −R , Fig. 2a. A ray emitted at an angle α is refracted at the interface. The angle of refraction β is determined by Snell's law: n 2 sin α = n 1 sin β . The refracted ray either crosses the x axis at some x c or disappears in the case of total internal reflection (TIR) if α > arcsin(n 1 /n 2 ) . The intersection can be virtual ( x c < 0 at β < 2α ) or real ( x c > 0 at β > 2α ). The regimes are determined by n r and α , see see Fig. 2c. In the limit α → 0 one can recover from (3) the paraxial approximation (2). The variation of x c (α) near the focus for paraxial rays x c (α → 0) is referred to as spherical aberration. For rather small n r , see n r = 1.3 in Fig. 2c, x c is negative and depends on α very weakly and, therefore, a sharp virtual image is expected to be formed. At n r = 1.8, 1.9, 1.95 , the virtual image position |x c | increases and its variation with α becomes very significant. This results in image blurring and, probably, disappearance. Moreover, at n r > √ 2 one obtains x c → −∞ as α → α * , with α * = arccos(n 2 /(2n 1 )) < π/4 . When α exceeds α * , the intersection x c becomes real and strongly dependent on α . At n r = 2.05, 2.1, 2.2 the dependence on α becomes flatter and that should increase the sharpness of the real image. At n r = 2 there is no intersection with x, i.e., x c → ∞ as α → 0.
(1) www.nature.com/scientificreports/ Thus, the rays emitted from the point source at x = −R do not converge to a well defined spot for 1.9 n r 2.2 . Under the severely strong spherical aberration the appearance of well resolved images is unlikely. Yet, the experiments show distinct images, see Fig. 1b. Thus, the formation of images and their locations observed experimentally seem to be in disagreement with the GO predictions.

Wave optics of ball lenses
Theoretical model. To explain the image formation let us now adopt the wave theory. We consider a 2-D model, in which the fields depend on x and y only, and take a point current source J z (x, y) = j 0 δ(x − R)δ(y) , see the coordinate system in Fig. 3b. Without the ball, the same point source located at x = y = 0 generates in the uniform medium an outgoing cylindrical wave with the electric field E u 0 is the Hankel function, ρ = x 2 + y 2 , k 1 = n 1 ω/c , and E 0 = j 0 ω/c 2 . Note that |E u z (k 1 ρ)| = E 0 at the distance ρ/ ≈ 0.9969/n 1 from the source so that E 0 can be conveniently used for normalization. The fields can be calculated without any further approximations by solving the Maxwell equations using the expansion into the cylindrical functions 23,26,27,38 . In general, modeling rigorously optical structures with sizes of hundreds and thousands of wavelengths thick is very challenging. Often this requires adopting various wave-propagation methods 39 which rely on some approximations, such as neglecting back-propagating fields. In the diffraction theory the field intensity is also represented as an integral over some aperture with a kernel which depends on the fields 40 . Since the fields at the aperture are not known, various assumptions are made, such as taking the incident field only (the Kirchhoff 's approximation) 41,42 . The specific circular geometry of the present problem allows finding the fields with accuracy limited only by computational precision rather than by a priori assumptions. The comparison of 2-D and 3-D FDTD results of modeling the focal distances of 10-18 µm-diameter spheres shows only a small difference while all trends remain essentially the same 43 . Thus, the application of a 2-D model here for significantly larger spheres is quite rational from the computational point of view and its results are expected to hold also in 3-D.
In the simulations we take a vacuum background n 1 = 1 and various values for n 2 . The results for n 1 = 1 can be obtained directly from those for n 1 = 1 as it follows from the Maxwell equations. Indeed, let us assume that we have some arbitrary current that produces fields in two situations, (a) and (b), so that the spatial dependences of the refractive index differ by a constant factor s: n b (r) = s n a (r) . The fields in case (b) are related to that in (a) as E b (r, ω) = E a (r, sω)/s and H b (r, ω) = H a (r, sω) . Thus, the case of a background index n 1 > 1 is equivalent to the case with n 1 = 1 but with operation at a proportionately higher frequency or, equivalently, shorter wavelength. This also increases the resolution by the same factor as in immersion microscopy.
Image plane position and magnification. The spatial spectrum of the fields in the far-field region from the ball at x ≫ R defines the image (real and virtual) which is formed by the objective (or cellphone camera) lens, see Fig. 1a. The image can be found by recreating the intensity in the focal plane of the objective lens using the far fields only 23,24,26,44,45 . This is equivalent to backpropagating the fields from some plane x ≫ R to the image plane. The distribution of the image intensity along y for a given location of the focal plane x gives the 2-D point spread function y) is the backpropagated field. The backpropagated field E b z (x, y) coincides with the true field E z (x, y) at distances larger than a few wavelengths from the ball at x > R . In the other regions they differ from one another. In particular, the true field diverges at the location of the point source while the backpropagated field (and PSF) does not. Note that without the ball, the z-polarized current point source located at the origin x = y = 0 produces far fields that after backpropagation give at x = 0 the image intensity psf(0, y) = π 2 E 2 0 J 2 0 (k 1 y) , where J 0 is the Bessel function. This implies the angular range of −π/2 < ϕ < π/2 for the collection of the outgoing waves. This image has the first zero at y/ = 2.4/(2πn 1 ) = 0.38/n 1 and its full width at half maximum (FWHM) is W/ = 0.36/n 1 . So the z-polar- www.nature.com/scientificreports/ ized source gives a smaller width compared to a y-polarized one 23 which can be attributed to its uniform angular emission.
To understand the role of wave phenomena in image formation we take a ball with R/ = 15 ( kR = 94.25 ) and n 2 = 2 , see Fig. 3. While its size is smaller than in the experiments shown in Fig. 1a, b it is still much larger than the wavelength. The paraxial GO predicts no ray convergence in this case, see Eq. (2). The intensity |E z (x, y)| 2 diverges near the point source and its value exceeds the color scale range in Fig. 3a. The intensity |E z (x, y)| 2 inside the ball near its boundary is also very high because of TIR. This accumulation of energy takes place even without any resonances, such as WGMs. The WGMs can also exist but their quality factor is extremely high around R/ = 15 with n 2 = 2 so that their excitation requires extremely fine tuning of R/ or n 2 . The psf(x, y) in Fig. 3b coincides with the true intensity |E z (x, y)| 2 in Fig. 3a for x > R but differs substantially for x < R . In particular, psf(x, y) has no enhancement inside the ball or singularity near the point source. The peak at x/ = 124 is the diffraction focus which is formed in contrast to the lack of ray convergence in the paraxial GO. Its location can be considered as the image plane x i . The longitudinal width of the peak (or the focusing depth) is very large, 20 . However, the transverse width is very small. Indeed, the FWHM of the central peak is W = 3.42 . Taking into account the magnification M = 124/15 = 8.27 , the resolution becomes W/M = 0.41 , which is close to the PSF width of 0.36 obtained in free space in 2-D for a z-polarized point source under the condition of light collection in the largest possible angular range of −π/2 < ϕ < π/2 . Thus, the ball itself does not produce any significant broadening of the PSF.
Having established the presence of focusing enabled by the contact ball with R/ = 15 , let us now move on to larger balls, closer to the GO regime. To investigate the location of the image plane we plot in Fig. 4 the PSF along the optical axis as the refractive index of the ball n 2 changes from 1.3 to 2.1 at a fixed size R/ = 100 . At n 2 = 1.3 a sharp peak is observed at x/R = −1.91 corresponding to the location of the virtual image, see Fig. 4a. As n 2 increases, the peak rapidly becomes smaller and broader. At n 2 = 1.9 , a noticeable intensity also appears at x/R > 0 . At n 2 = 1.95 the furthest maximum of intensity is at x/R ≈ 50 , however, the second one at x/R ≈ 15  www.nature.com/scientificreports/ is higher, see Fig. 4b. As n 2 increases even more, the furthest maximum becomes narrower, higher, and closer to the ball. The furthest peak is still slightly smaller than the second peak at n 2 = 1.98 but it overcomes it at n 2 = 2 . Note that psf(x, 0) always oscillates along x. These oscillations can disappear in practice if imaging is performed using illumination with a broad spectrum. However, the slow variation can still remain if the illumination bandwidth is sufficiently small. The plane locations at the maxima along x that differ from the furthest one cannot produce images because the PSF has much stronger sidelobes in the transverse y direction, see, for example, the plane at the second maximum at x/ = 56 in Fig. 3b. Thus, the furthest peak defines the image plane location x i . We can now compare the predictions of GO with that of WO. Due to symmetry of the problem, M in both cases is directly related to the image location M = −x i /R . Since x i and M diverge as n r → 2 in the GO limit, it is convenient to analyze graphically R/x i = −1/M . Figure 5a,b show the comparison for virtual and real imaging, respectively, for several ball sizes: R/ =50, 100, and 1000. For virtual images, WO predicts slightly larger |x i |/R (or smaller R/|x i | ) as compared to that from GO. However the agreement with GO becomes closer as R/ increases. More interesting situation is observed for real images. The function R/x i for all R/ behaves linearly with n 2 having practically the same slope as predicted by GO. However, the lines are shifted and the shift depends on R/ . With increasing R/ , the shift decreases and eventually one recovers the GO regime. The presence of even a small shift in R/x i in the region of n r ≈ 2 (or 1/M ≈ 0 ) can lead to dramatic deviations of x i /R and M from the GO predictions. The image plane positions calculated using WO for the indices n 2 and sizes R/ realized in the experiments are given in Table 1 and agree well with the experimental measurements. Figure 5c shows magnification (and therefore, the position of the real image plane x i = |M|R ) as a function of R/ for n 2 = 2.02 , 2.05, and 2.1. In all cases |M| initially grows with R/ . If n 2 differs from 2 significantly, for example, n 2 = 2.1 , then |M| saturates quite rapidly at R/ ∼ 200 − 300 . If the difference n 2 − 2 becomes smaller, then the growth of |M| becomes much slower and reaching the asymptotic GO values requires much larger values of R/ . For n 2 = 2.05 , for example, |M| still continues growing slightly even for R/ ∼ 1000 . For n 2 = 2.02 , the growth of |M| is very pronounced even for R/ ∼ 2000 . Thus, the convergence to the GO regime is determined not only by the ratio R/ but also by the difference n 2 − 2 . When the difference n 2 − 2 is small, the deviation from the paraxial GO even for large lenses remains very significant. This conclusion is supported well by the experimental results, see Table 1.
Resolution. According to the Houston criterion the resolution is defined as the FWHM of the PSF 46 . After finding the image of the point source, the width W of the central peak was divided by |M| and normalized to . Figure 6a shows the transverse dependence of psf(x i , y) at the image plane x i for various n 2 and fixed R/ = 100 . In all cases the PSF has sidelobes. Their height relative to the central peaks decreases as n 2 increases from 1.98 to 2.05. In the estimates of the FWHM of the PSF we only consider the central peak keeping in mind the adverse effect of large sidelobes. The FWHM of the central peak does not change significantly: W/|M| = 0.64 at n 2 = 1.98 and W/|M| = 0.67 at n 2 = 2.05 . Although this resolution slightly worse than the diffraction limit /2 , it is surprisingly high considering the large spherical aberration, see Fig. 2c. Thus, for a fixed R/ the resolution does not change significantly for different n 2 but magnification at n 2 = 1.98 is larger than that at n 2 = 2.05 . On the other hand, the larger sidelobes at n 2 = 1.98 will lead to lower image quality. The resolution decreases with increasing R/ . For example, in Fig. 6b the FWHM increases from W/(|M| ) = 0.98 at R/ = 500 to W/(|M| ) = 1.38 at R/ = 2000. Figure 6c shows that the PSF width W/(|M| ) for a fixed n 2 grows with increasing R/ but does not significantly depend on n 2 . On the other hand, the magnification |M| increases both with n 2 and R/ , see Fig. 5c. Thus, for a fixed n 2 the increase of |M| always comes at the expense of decreasing resolution. This trend is illustrated in Fig. 6d which shows the relation between the resolution and magnification for several n 2 . A natural question arises of whether we can achieve simultaneously high resolution (small PSF width) and large magnification.
Unfortunately, the answer to this question is negative. Indeed, let us fix a rather high |M| = 40 and try to increase where r ′ is a point at the cylindrical surface with radius R. In Eq. (4), the effective electric and magnetic currents, which are obtained from the fields outside the cylinder 47 , and the Green's functions are where H 0,1 are the Hankel functions. To investigate the contribution to E z (r) from different parts of the integration surface, let us define the cumulative field, which is formed by the currents in the angular range limited by [−ϕ : ϕ], so that the integral over the full circle ϕ = π gives the actual field at the observation point: E z (r) = E c z (r, π) . Note that the field is calculated using the integration over a closed surface, not a finite opening in a screen which is common in the diffraction theory. Figure 7a shows the magnitude of the cumulative field |E c z (x, ϕ)| as a function of ϕ for n 2 = 2.02 at several special locations on the x axis: x/R = 101 (GO focus), x/R = 20.3 (diffraction focus), and x/R = 12.77 (the furthest minimum of PSF). At the GO focus, the contributions |E c z (r, ϕ)| originate from a rather small angular range ϕ/π 0.15 . At the diffraction focus, the contributions increase more rapidly, come from a larger interval ϕ/π 0.25 , and subsequently result in a much larger value for |E z (r)| . At the PSF minimum, the contributions initially increase and then start to decrease leading eventually to a very small |E z (r)|.
The rigorous representation of the field at any observation point using the fields at the circle just outside the ball, see Eq. (4), can also be recast into a simplified form. We can assume that the effective currents are The term f 2 describes defocus aberration and f 4 describes spherical aberration 48 . By setting f 2 = 0 in Eq. (8) we can obtain the focal distance d f = 2R/(n r − 2) and magnification in the paraxial GO approximation (2). The focusing properties of an optical systems are often characterized by the Fresnel number N of the exit pupil with radius R: N = R 2 /( d f ) , where d f is the GO focal distance. Typically, the presence of a sharp focus requires N ≫ 10 . For N 10 , the defocus tolerance can be significant. Moreover, the maximum of irradiance can shift closer to the exit pupil. The regime N 10 can significantly differ from the predictions of GO due to the influence of diffraction which defines not only the resolution but also the position of maximum intensity. In the case of ball lenses, there is no clearly defined aperture. However, to estimate the Fresnel number we can simply use the ball radius R since the effective aperture is formed by the TIR: For large d f , which are obtained for small n r − 2 , the Fresnel number becomes small. For example, the Fresnel number reaches N = 10 only at R/ = 1000 if n r = 2.02 . As seen from Fig. 5c, at R/ = 1000 the focal position is still significantly smaller than the GO prediction.
The phase difference allows us to explain the shift of the intensity maximum as compared to the paraxial GO. Figure 7b shows the phase difference θ(α) as a function of the location on the circle ϕ = 2α for the same observation points on the x axis as in Fig. 7a. In general, if the spherical aberration term vanishes f 4 = 0 , the focal point is defined by f 2 = 0 and its intensity is determined by the integral over the full aperture (limited by TIR) resulting in the largest possible intensity. The deviation f 2 = 0 from the GO focus leads to the blurring of the image due to defocus aberration. The spherical aberration f 4 = 0 , see the line for x/R = 101 in Fig. 7b, gives rise to a rapid phase decrease for α/π > 0.15 and therefore, to the oscillations of the integrand reducing significantly its value. A larger intensity can be obtained if one moves away from the GO predicted focus in such a way as to provide the smallest variation of the phase θ(α) over the aperture. Since usually f 4 < 0 , compensating it requires f 2 > 0 , and therefore, taking smaller d as compared to the paraxial GO prediction. Indeed, at the point of maximum, x/R = 20.3 in Fig. 7b, the phase deviation from zero is limited to |θ| π/2 in the interval ϕ/π < 0.25 . For even smaller x, the minimum at x/R = 12.77 in Fig. 7b, the phase changes significantly bringing the total intensity to almost zero.
We can further apply this simple physical picture to derive an analytical curve that describes the image plane location and magnification. From Fig. 7b we can note that at the location of intensity maximum the phase www.nature.com/scientificreports/ difference is limited to θ π/2 . Indeed, exceeding this value would cancel out the contributions from the smaller angles. Thus, we can state that at the extremum θ(α * ) = π/2 . Using (8) we can find α * and then relation between f 2 and f 4 . Further, we can approximate f 4 ≈ −1/2 and obtain The paraxial GO result described by Eq. (2) follows from the more general Eq. (10) in the limit R/ → ∞ . The WO effects in the position of the image plane and corresponding magnification are accounted for by the extra term δ , which depends on the size parameter kR. Note that one can often use GO formulas and account for the effects not described by GO by introducing some effective parameters, for example, an effective refractive index of a spherical lens which becomes a function of its diameter 43 . In our case, we obtained the GO result from the more general WO result and, therefore, we do not need to resort to any effective parameters. Figure 5c shows that simple formula (10) describes the fully numerical results very accurately. This extra term also agrees with the kR-dependent shifts of the straight lines in Fig. 5b, especially for large R/ . Moreover, Eq. (10) explains accurately the experimentally measured image plane locations, see Table 1. The analytical derivation of Eq. (10) uses the phase difference between the path characterized by α = 0 and that for α = 0 . We neglected the variation of intensity as a function of α in order to obtain a simple result which, however, agrees well with the fully numerical modeling. In the 3-D case, the formula for the optical path difference remains valid but one may need to account for the solid angle factor sin α and specific orientation of the emitter.
We can conclude that there are two main physical factors that contribute to the shift of the image plane in contact-ball imaging. First, even an aberration-free converging spherical front does not focus to its origin once an aperture is introduced 41,42 . The focal point moves towards the aperture and the shift increases with decrease of the Fresnel number of the aperture. When the relative refractive index of the ball is n r ≈ 2 , the geometric focus goes to infinity and, therefore, the Fresnel number decreases drastically. Second, in the case of contact ball lenses the phase front immediately after the ball suffers from spherical aberration and, therefore, it cannot focus to a single point. In terms of ray optics, the spherical aberration is particularly significant near n r ≈ 2 when no single focal point for the rays refracted by the ball can be defined. Both factors define the distribution of the diffracted intensity in contact-ball imaging. Within the realm of WO the intensity at a given point can be written as a diffraction integral over an aperture effectively defined by TIR of the spherically aberrated wavefront. The axial intensity reaches a maximum at the point at which the defocus and spherical aberration terms, see Eq. (8), can partially compensate each other over the integration aperture. www.nature.com/scientificreports/ Experimental resolution. The resolution of the ball lens imaging can also be estimated from the experiments shown in Fig. 1a,b. The image of the Siemens star was scanned at several distances from the center, see Fig. 8a,b. The metallic spokes manifest themselves as the dips. The oscillations of the intensity become larger as the local period of the structure increases. The image intensity contains some spatially uniform background level I b due to various scattering in the experimental setup. Using the SEM of the Siemens star, see Fig. 1b, the gap-topitch ratio L 1 /L can be estimated as 0.6 in the area of the scans.
To extract the resolution from the measurements it was assumed that at each scan the star can be described as a 1-D structure in which each period L contains a bright stripe with intensity I 0 and width L 1 . The image is obtained by the convolution of this ideal object with a Gaussians PSF defined by its FWHM W. According to the Houston resolution criterion 46 , the FWHM of the PSF represents the resolution of the system. The resultant image is also a periodic structure which can be characterized by the intensity at the peak I peak and at the dip I dip . For a specific geometrical parameter L 1 /L , the intensity depends only on the ratio W/L and therefore, one can define where 0 < f (W/L) < 1 . The difference I peak − I dip does not depend on the presence of the background level I b in the measurements and, therefore, can be directly used for fitting once the geometric function f(W/L) is known. Figure 8c shows the function f(W/L) at several values of L 1 /L . Note that this function is symmetric relative to L 1 /L = 0.5 , that is it takes the same value, for example, for L 1 /L = 0.4 and L 1 /L = 0.6 . In all cases the modulation amplitude f(W/L) decreases with increasing W/L. The oscillations become observable for W/L 1.
The fitting of the analytical dependence (11) with L 1 /L = 0.6 to the experimental points, see Fig. 8d, gives W = 0.92 µ m or W/ = 1.9 . Varying W near this optimal value makes the fit worse and therefore, the approach allows a reliable extraction of resolution. Indeed, the curve for W = 1.02 µ m cannot be fitted to the points better than shown since any variation of I 0 would only scale the curve up or down and therefore, worsen the fit. The value W/ = 1.9 is slightly greater than the PSF obtained in the simulations W/(|M| )=1.4, see the curve for n 2 = 2.05 at R/ ≈ 2000 in Fig. 6c. The knowledge of I 0 allows estimating the background level I b as well. Without the background scattering and for small periods the intensity should oscillate symmetrically around its average value (L 1 /L)I 0 = 0.6I 0 = 0.55 while the average value for scan 1 is I ave = 1.19 . This gives the background level I b = I ave − 0.6I 0 = 0.64 . Thus, the background scattering is very significant and should be accounted for in resolution estimations.
We also note that for a gold double-stripe object the resolution was previously estimated 31 to be W = 0.65 µ m at = 480 nm or W/ = 1.4 . This resolution is better than in the present study of the star. The difference can be attributed to several factors. One factor is related to geometry since the double-stripe element of the long-period array in Ref. 31 is practically an isolated object and, therefore, is more suitable for resolution quantification 49 compared to the star with azimuthally periodic spokes, which also have a width gradient in the radial direction. Another factor is related to different reflection properties of thin metallic layers, Au in Ref. 31 versus Cr in the present study.

Conclusion
Despite a tremendous interest in MSI methods developed with wavelength-scale microspheres 5-21 , the connection of this field of studies to the standard wave theory of aberrated imaging by macroscopic ball lenses has not been previously established. In this work we developed a comprehensive approach to this problem for ball lenses with diameters varying from D/ ≈ 30 (quite often used in MSI studies) up to D/ ≈ 4000 (reaching the limit of millimeter-scale ball lenses). Our approach is based on the transition from geometrical optics to full-wave solutions of the Maxwell equations. A unique feature of our numerical modeling approach is that it bridges wave phenomena taking place at completely different spatial scales -the near-field coupling of a point source, the field propagation inside ball lenses with diameters 30 D/ 4000 , and, finally, the formation of the diffracted field at distances ∼ 10 5 .
Our theory is developed specifically for the case of a direct contact of a ball lens with nanoscale objects with intention to increase the resolution due to near-field effects in a spirit of the solid immersion lens concept 50,51 . Another feature of the theory is its ability to describe accurately the imaging by ball lenses with refractive index near the critical value n ≈ 2 , for which the deviations from geometrical optics become dramatic. This critical regime is very attractive for imaging applications because of extremely large M > 50 magnification 31 . We showed that the image plane location and magnification in this regime can be described correctly only by wave optics which predicts a significant shift of the image plane towards the ball with corresponding reduction in magnification, in accord with the presented experimental evidence. The shift becomes more significant as ball's size becomes smaller and n closer to 2. This effect is governed mostly by diffraction since the Fresnel number of the effective aperture is small. It is demonstrated that the exact location of the image plane is defined by the counterplay of defocusing and spherical aberrations. This allowed us to derive a simple correction to the geometric formula for the image plane location. In contrast to common expectation, we showed that the resolution improves for smaller ball lenses reaching deeply subwavelength values ∼ /2 at R/ 100 in the vicinity of n ≈ 2 . The index variation about n ≈ 2 affects greatly the magnification but not the resolution.
Although the dispersion n( ) is present in all microspheres used previously for MSI, its impact on imaging becomes dramatic only in the vicinity of the critical value n ≈ 2 for contact microspheres in the air environment studied here. These effects are significantly less pronounced for the relative indices 1.4 < n < 1.8 which are typical for the previously studied silica microspheres in air or BTG microspheres in water, in PDMS, or in plastic environments.
To experimentally test our theoretical predictions, we performed a resolution quantitation aimed at the demonstration of potential advantages of ball lenses with the specially designed index for cellphone-based www.nature.com/scientificreports/ microscopy applications. Currently, the resolution of such portable and lightweight microscopes is pixel-limited at a ∼ 1.5-µ m level due to insufficient magnification. Using a Siemens star object imaged through a LASFN35 glass ball lens with n = 2.049 at = 480 nm we demonstrated a resolution of ∼ 0.9 µ m or 1.9 , slightly below the resolution of 1.4 predicted by our theory, and magnification M ≈ 35 . In our previous studies we approached the wavelength-scale resolution using similar approaches 31 . Thus, this experiment demonstrates a potential of the proposed methods for increasing the resolution of imaging based on using contact ball lenses with n ≈ 2 . The peculiar imaging properties of contact ball lenses with n ≈ 2 complement their ability to focus plane waves to strongly enhanced photonic nanojets on the outer edge 52 . The theoretical methods developed in this work demonstrate a gradual transition from MSI methods to classical imaging by aberrated macroscopic ball lenses with applications in high-resolution cellphone-based microscopy.

Methods
Experiments. All experiments were performed using the cellphone-based transmitted light microscopy shown in Fig. 1a. The illumination from a tungsten halogen lamp was provided through narrow (about 10 nm) bandpass filters with transmission peaked at different wavelengths. A white light 120-grit ground glass diffuser (Edmund Optics) was installed at 8-mm distance below the sample to provide widefield incoherent illumination with a broad range of incident angles. The object was a Siemens star (Edmund Optics) made of 36 equidistant Cr spokes on a silica substrate. The ball lens was made of LASFN35 glass with the refractive index varying from about 2 to 2.1 as the wavelength changes from 700 nm to 400 nm. The measurements were performed in air environment. The optical images in Figs. 1b and Fig. 8a were recorded using a cell phone camera (Samsung Galaxy S9+). The SEM image in Fig. 1b was obtained using the Raith 150 e-beam lithography system. The location of the scans in Fig. 8a were selected close to the central part of the image to minimize pincushion distortion. The center of the Siemens star was shifted to the edge of the limited circular field-of-view where all distances are locally distorted. The local coordinate scale along lines (1-3) in Fig. 8a was determined by the magnification data independently obtained for the same ball lens at the central part of the image using a double-stripe object with known physical dimensions.
Simulations. All simulations were performed using in-house developed codes in C programming language.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.